Normal_t.yb <- function(data, Col_DV, Col_IV, Col_group = NULL) {
  library(psych)
  library(plyr)
  library(car)
  df <- data
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  df <- unite(df, "IV", all_of(Col_IV), sep = "_", remove = F)
  df$IV <- paste("IV", df$IV, sep = "_")
  df$IV <- as.factor(df$IV)
  out <- list()
  out[[1]] <- df
  out[[2]] <- data.frame(group = "", DV = "", Conversion_process = "")[-1,]
  for (i in c(1:length(group))) {
    df1 <- df[df$group == group[i],]
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    for (j in c(1:length(Col_DV))) {
      homogeneity <- c()
      for (r in 1:4) {
        if (r == 1)
          df1$y <- df1[,Col_DV[j]]
        if (r == 2)
          df1$y <- df1[,Col_DV[j]]^0.5
        if (r == 3)
          df1$y <- log(df1[,Col_DV[j]])
        if (r == 4)
          df1$y <- 1/df1[,Col_DV[j]]
        homogeneity[r] <- as.data.frame(leveneTest(eval(parse(text = paste("y", paste(Col_IV, collapse = "*"), sep = "~"))), data = df1, center = mean))[1,3]
      } 
      if (homogeneity[1] > 0.05)
        Homogeneity_variance <- "Homogeneity"
      else {
        if (any(homogeneity[2:4] > 0.05)) {
          r <- which(homogeneity[2:4] == max(homogeneity[2:4]))+1
          if (r == 2)  {
            df1[,Col_DV[j]] <- df1[,Col_DV[j]]^0.5
            Homogeneity_variance <- "+SQRT"
          }
          if (r == 3) {
            df1[,Col_DV[j]] <- log(df1[,Col_DV[j]])
            Homogeneity_variance <- "+Ln"
          }
          if (r == 4) {
            df1[,Col_DV[j]] <- 1/df1[,Col_DV[j]]
            Homogeneity_variance <- "+1/x"
          }
        } else {
          Homogeneity_variance <- "Non homogeneity"
        }
      }
      Is_normal <- F
      r <- 0
      while (all((!Is_normal), r < 10)) {
        r <- r + 1
        normal <- data.frame(group = "", DV = "", IV = "", n = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
        for (k in 1:length(IV)) {
          sample <- df1[df1$IV == IV[k], Col_DV[j]]
          n <- length(sample)
          ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
          sw <- shapiro.test(sample)
          normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), n, 
                                    round(ks[[1]], 4), round(ks[[2]], 4), 
                                    round(sw[[1]], 4), round(sw[[2]], 4)))
        }
        names(normal) <- c("group", "DV", "IV", "n", "KS", "KS_P", "SW", "SW_P")
        if (mean(as.numeric(normal$n)) < 100) {
          Is_normal <- all(normal$SW_P > 0.05)
        } else {
          Is_normal <- all(normal$KS_P > 0.05)
        }
        if (!Is_normal) {
          Skewness <- data.table::rbindlist(
            describe(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), idcol="IV")
          Skewness$skew.se <- sqrt(6 * Skewness$n * (Skewness$n - 1)/(Skewness$n - 2)/(Skewness$n + 1)/(Skewness$n + 3))
          Skewness$skew.Z <- Skewness$skew/Skewness$skew.se
          if (mean(Skewness$n) <= 50) 
            z_0.05 <- 1.96
          if (mean(Skewness$n) <= 300 & mean(Skewness$n) > 50) 
            z_0.05 <- 3.29
          if (mean(Skewness$n) > 300) 
            z_0.05 <- mean(2/Skewness$skew.se)
          if (abs(mean(Skewness$skew.Z)) < z_0.05) {
            if (r == 1)
              Conversion_process <- "skew ∈ normal"
            else
              Conversion_process <- paste(Conversion_process, " skew∈normal", sep = " ~")
            Is_normal <- T
          } else {
            Conversion_process <-""
            if (mean(Skewness$skew) > 0) {
              if (abs(mean(Skewness$skew.Z)) < 3*z_0.05)  {
                df1[,Col_DV[j]] <- df1[,Col_DV[j]]^0.5
                Conversion_process <- paste(Conversion_process, " +SQRT", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 3*z_0.05 & abs(mean(Skewness$skew.Z)) < 7*z_0.05) {
                df1[,Col_DV[j]] <- log(df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " +Ln", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 7*z_0.05) {
                df1[,Col_DV[j]] <- 1/df1[,Col_DV[j]]
                Conversion_process <- paste(Conversion_process, " +1/x", sep = " ~")
              }
            }
            if (mean(Skewness$skew) < 0) {
              if (abs(mean(Skewness$skew.Z)) < 3*z_0.05)  {
                df1[,Col_DV[j]] <- (max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])^0.5
                Conversion_process <- paste(Conversion_process, " -SQRT", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 3*z_0.05 & abs(mean(Skewness$skew.Z)) < 7*z_0.05) {
                df1[,Col_DV[j]] <- log(max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " -Ln", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 7*z_0.05) {
                df1[,Col_DV[j]] <- 1/(max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " -1/x", sep = " ~")
              }
            }
          }
        } else {
          if (r == 1)
            Conversion_process <- "Normal"
          else
            Conversion_process <- paste(Conversion_process, " Normal", sep = " ~")
        }
      }
      if (r == 10)
        Conversion_process <- paste(Conversion_process, " 10times", sep = " ~")
      out[[1]][(out[[1]]$IV %in% df1$IV & out[[1]]$group == df1$group), Col_DV[j]] <- df1[,Col_DV[j]]
      out[[2]] <- rbind(out[[2]], c(group[i], Col_DV[j], Homogeneity_variance, Conversion_process))
    }
  }
  names(out[[2]]) <- c("group", "DV", "Homogeneity_variance", "Conversion_process")
  out[[1]] <- out[[1]][,colnames(data)]
  return(out)
}  #正态性检验；输入(列名)：data-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
   #输入：out[[1]]：转换后的数据与输入的data相同；out[[2]]：转换过程(Homogeneity_variance=方差齐性；Conversion_process=正态性转换)

   #示例
data <- data.frame(IV_A = as.factor(rep(c(rep(1,20),rep(2,20),rep(3,20),rep(4,20)),2)),
                   IV_B = as.factor(rep(rep(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5)),4),2)),
                   re = rep(rep(c(1:5),16)),
                   group = as.factor(c(rep(1,80),rep(2,80))),
                   DV_1 = runif(160),
                   DV_2 = rnorm(160))
                           #创建示例数据框，IV_A&IV_B(自变量A、B, 需为因子factor), re(重复), group(分组), DV_1&DV_2(因变量1、2)
data                       #查看示例数据框

out <- Normal_t.yb(data, c("DV_1", "DV_2"), c("IV_A", "IV_B"), c("group"))